In order to write a computer program to calculate this continuous approximation of the terrain, we need to do the following.
First, we have to describe our problem in a formal, mathematical way. In particular, we have to choose the mathematical definition of the object - bitmap, and the mathematical definition of our continuous description of the terrain represented by the bitmap.
It seems reasonable to say that a bitmap is a function defined in an area
\( \Omega = [1,maxx]\times[1,maxy] \ni (x,y) \rightarrow BITMAP(x,y) \in [0,255] \)
where by \( \Omega \) we mark the entire area where our bitmap is spread.
In turn, our continuous representation of the world will be represented by a function
\( u \) of real values.
\( \Omega = [1,maxx]\times[1,maxy] \ni (x,y) \rightarrow u(x,y) \in [0,255] \)
We want our function to be "smooth" and continuous. Mathematically, we write the condition that our function will be of the
\( C^1 \) class, it means that at each point it will be smooth enough to compute its derivatives in directions perpendicular to the edges of our area. This practically means that at every point we can apply a "ruler" to the graph of our function perpendicularly to one of the edges of our area, and measure the angle between this ruler and the base (a flat surface spread out at zero height). After all, the derivative is nothing but the tangent of this angle. In other words, our function will have no "kinks" or gaps where you would not know how to use our ruler. At the bends, this ruler would jump from one position to another, and in the case of a hole, it would not be clear how to use it. Of course, the ability to measure the derivative (placing a ruler) in two directions perpendicular to the edge of the area also means the ability to apply a ruler and measure directional derivatives in any other directions not perpendicular to the edge. Hence, we can smoothly move across the plot of such a smooth function in any direction.
How to get such a continuous smooth function? We have to decide how our function will be constructed. For example, we can divide the area over which the bitmap spans into certain elements, and based on these elements we can define a set of many smooth functions, from which we then "glue" our \( u \) function.
Let us imagine that at the height corresponding to the zero height (equal to the value of the zero pixel) Let us imagine that at the height corresponding to the zero height (equal to the value of the zero pixel) we build a flat two-dimensional mesh, of which the number of square elements is arbitrary. These meshes, according to the adopted nomenclature, will be called finite elements, because each of these elements has a limited finite area. We might have fewer elements than the pixels in the bitmap, and then several pixels will be spread over each such element. The boundaries between our elements do not have to match the boundaries of the pixels. They can be freely defined on a flat surface. It may also be the case that there are more elements than pixels, and then there will be many such elements in each pixel. Let us assume, however, that we have fewer elements than pixels. They form a regular grid of \( N_x*N_y \) finite elements.
Now, on each such element, we define a smooth function. We can use B-spline functions for this purpose.
These functions were first introduced by the American mathematician of Romanian origin, Isaac Jakub Schoenberg [1].
B-spline functions are widely used in computer modeling and simulations thanks to the growing popularity of the field called isogeometric analysis disseminated by prof. T. J. R. Hughes [2].
The idea behind these methods is to use families of B-spline functions for finite element calculations.
We denote these functions by \( B_{i,j;2}(x,y) \), where \( i \) and \( j \) denote the numbering of our functions, and 2 means that they are polynomials with pieces of the second order, class \( C^1 \).

Fig. 1 shows three exemplary two-dimensional B-spline functions spanned over a two-dimensional mesh. Each such two-dimensional B-spline function is created by selecting and multiplying by itself two one-dimensional B-spline functions, one selected from the set of one-dimensional B-spline functions spanning along the horizontal axis of the mesh, and the other selected from the set of one-dimensional B-spline functions spanning along the vertical axis of the mesh. These sets are called the basis of one-dimensional B-spline functions.
In turn, we denote these one-dimensional B-spline functions \( B^x_{i;2}(x) \) and \( B^y_{j;2}(y) \) where variables \( x \) and \( y \) identify the directions (coordinate axes) along which our B-splines (horizontal axis \( x \) and vertical axis \( y \)), \( i \) and \( j \) numbering areas of these functions (which next one-dimensional B-spline side we choose from such one-dimensional basis), and 2 means again that they are polynomials of the second degree, of the class \( C^1 \) (i.e. that we can compute the first derivatives of them).
Then the selected one-dimensional functions are multiplied by themselves, which gives us a smooth two-dimensional B-spline basis. Such a relaxation of a function by multiplying two-dimensional functions is called a tensor product of a function.
\( B_{i,j;2}(x,y)=B^x_{i;2}(x)B^y_{j;2}(y) \)
This is illustrated in Fig. 1 and Fig. 2.
The two-dimensional B-spline functions given in this way have the shape of a slender "hill", with similarity among nine adjacent elements. The highest point of such a function - the hill - is in the center of the middle element. These functions go down to the zero value, assumed at the edges of the span defined by the nine elements on which the function is defined. These functions, according to the nomenclature, are called basis functions.

Our continuous approximation of the terrain is obtained by adding together many such slender hills - B-splines. Each of them will be scaled (raised up or down) so as to obtain a total of the continuous approximation of the terrain. If we correctly select the heights of individual hills, then we will get smooth approximations of our terrain as presented in Exemplary problem of two-dimensional bitmap projection-Fig. 2.
Now the question is how to choose the heights to which we will draw our basis functions. The first method that comes to mind is to choose the pixel value z \( BITMAP(x,y) \) exactly at the highest point of the B-spline (in the middle of the hill). Unfortunately, this method has several disadvantages. First, if there are fewer B-spline basis functions than pixels, then we ignore all neighboring pixels in our B-spline area, choosing only one value from the center of the area. Our bitmap is likely to have some disturbances and the fact that we will hit a localkickback is a measurement error, or we may hit a local hole in the terrain, or a local tree or building that has disturbed the measurement of the terrain. Secondly, we will notice that our basis functions extend over a square of nine elements. Since each such B-spline function can be associated with the center of its element, and each of these functions extends over nine adjacent elements, it means that there is a total of nine such functions spread on each element, and that the adjacent functions overlap. So, if we were to stretch the functions so that their maximum point coincides with the central pixel, and we summed all these functions together to get our global function \( u \), then on each element, even in the central point, our resulting approximation (the sum of these functions) would be higher than our central pixel. It is because the adjacent eight functions would also be non-zero on a given element and would raise the value of our approximation at this point up.